Conversation
…-limit bypass for testing Wind coupling (opt-in, backward compatible): The default directional projection formula evaluates Rothermel's model with the wind projected onto the front normal before raising it to the B power, R(theta) ~ (U cos theta)^B. This is non-convex once phi_w(B-1) > 1, so a finite ignition line's endpoints grow a Wulff wedge whose tip falls well below the head rate over time -- verified against a Munoz-Esparza-style finite-line test case, where the head rate degraded from ~100% to ~40% of the theoretical Rf over 2100s. erf.fire.directional_wind_coupling = "wrf" (default: "projection", unchanged) instead exponentiates Rothermel's wind/slope factor from the raw, unprojected wind speed and scales the result by the cosine to the spread direction afterward, matching WRF-Fire's fire_ros (module_fr_fire_phys.F). This keeps R(theta) linear in cos(theta) -- the convex support function of a circle -- so the wedge does not form; the same test case now tracks the theoretical head rate to within 0.5% for the full run. Default behavior is verified bit-for-bit identical to the prior code (max|delta phi| = 0 over the whole domain on the same test case). Test infrastructure: - erf.fire.prescribed_wind[_x/_y]: bypass atmospheric wind interpolation with an exact constant wind, for controlled idealized-wind test cases. - erf.fire.use_wind_limit: allow disabling the MEWS wind-speed cap (default on, unchanged), needed when a test's wind speed exceeds the cap and the cap would otherwise mask the ROS model's true behavior. Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com>
hgopalan
left a comment
hgopalan
left a comment
There was a problem hiding this comment.
Thanks for this. I audited it in detail: test-merged onto current ERF-Hazard (clean, no conflicts), built Release with unit tests on (zero warnings, all 20 gtests in the DirectionalShape / WindLimit / LevelSetAdvection suites pass), and ran the existing Exec/RegTests/FireDirectionalShape decks with base and merged binaries.
Verdict: the WRF coupling itself is correct and does what the description says. But the PR also bundles a second, undocumented change to the front-normal stencil that alters every default directional run, so the "bit-for-bit identical" claim does not hold for the default scheme. Requesting changes to split that out.
Blocking
1. Hidden change to the front normal
fill_directional_ros now computes the normal as 0.5*(dm + dp) of one_sided_derivatives with the run's gradient scheme, instead of central differences. The default scheme is weno5z_front, so inside the front band the normal comes from two HJ-WENO5-Z reconstructions: a different stencil. This is not mentioned in the title, description or commit message.
Existing FireDirectionalShape projection deck, base vs merged binary, 1500 s:
| deck | base head [m/s] | merged head [m/s] | arrival times |
|---|---|---|---|
projection, default weno5z_front |
0.18487 (40 % of the way from the Wulff tip to Rothermel) | 0.14949 (7 %) | differ everywhere |
projection, erf.fire.levelset.gradient = upwind |
0.15279 | 0.15279 | identical (phi differs by 1e-9, rounding of the averaged stencil) |
So the bit-for-bit claim holds only on an upwind deck; the test that produced "max|delta phi| = 0" must have used levelset.gradient = upwind. The sharper normal pushes the projection closer to its exact Wulff tip, which is defensible in isolation, but it changes the README table of that regtest and every default directional run. It also leaves the code inconsistent: the hybrid path (ERF_HybridRos.H) calls fill_directional_ros positionally and so keeps a first-order normal, and Balbi keeps central differences.
Please revert the normal change in this PR and open it separately, with the re-measured FireDirectionalShape README table and the three places that still say "central differences" (file comment in ERF_DirectionalRos.H, Docs/sphinx_doc/theory/fire_propagation.rst, the checker docstring) updated.
2. Stale description
The description and commit message list erf.fire.use_wind_limit as new test infrastructure. That flag is already in ERF-Hazard from #381 and this diff contains no change to it.
Should fix
3. Validate the option combination at parse time
directional_wind_coupling = "wrf" is silently ignored with directional_shape = "ellipse" and with every model other than Rothermel. The surrounding parse block aborts on contradictory combinations (Anderson flanks without the ellipse), so this should abort too. BEHAVE (behave_ros_cell) uses the identical power-law form, so extending the coupling to it is a few lines and removes one of the silent cases.
4. No regression test
The finite-line case that motivated the change is not in the PR, and the convention here is opt-in default plus a deck that compares old and new. Cheapest route: add inputs_wrf (the projection deck plus erf.fire.directional_wind_coupling = "wrf") to Exec/RegTests/FireDirectionalShape and a wrf kind to check_firedirectionalshape.py whose exact normal speed is R0 (1 + phi_w(U) max(cos theta, 0)). I did that in a scratch copy with the merged binary:
wrf:
head (+x) PASS 0.24846 m/s from 181 cells vs Hopf 0.24945 m/s (-0.39 %)
head vs Rothermel PASS 0.24846 vs 0.24943 m/s (-0.39 %)
back (-x) PASS 0.02396 m/s vs Hopf 0.02405 m/s
flank (+y)/(-y) PASS 0.02403 m/s vs Hopf 0.02403 m/s
That is a one-deck proof of the claim on a point fire. The finite-line deck would be a welcome second test, since FireLineFire is periodic in y and so an infinite line.
5. Documentation and decks
No rows for the four new keys in Docs/sphinx_doc/Inputs.rst, nothing next to the directional_shape paragraph in the theory doc, and Exec/CanonicalTests/Fire/inputs_fire_master_reference does not list them. rothermel_ros_cell also needs @param entries for the two new cosine arguments.
6. Position the option against the existing fix
directional_shape = "ellipse" already removes the wedge and gives -0.67 % on the same deck. The file comment presents the ellipse as the only remedy; it should say what "wrf" adds (the model's own oblique rates, linear in cos theta, rather than an imposed ellipse, and the form WRF-Fire users will compare against).
Minor
- The comments call
R0 (1 + phi_w cos theta)"the support function of a circle". With the clamp at zero for backing it is the support function of a stadium (a disc of radius R0 swept along the wind vector). Still convex, so the conclusion stands, but the wording is off. prescribed_windbypasses the interpolation cleanly andfire_wind_refhas no ghost cells, so thesetValis fine. WAF and the FARSITE terrain corrections still apply downstream, so "exact constant wind" is only true withuse_waf = falseanduse_terrain_wind = false; worth saying in the parameter comment.- Checked against WRF-SFIRE
fire_ros(module_fr_fire_phys.F): it caps the raw speed at 30 m/s, exponentiates, then multiplies bymax(0, (v.n)/|v|)andmax(0, (dz.n)/tanphi). The PR matches this apart from the 30 m/s cap, which ERF covers with the MEWS limit.
🤖 Generated with Claude Code
Runs a finite ignition line under a strong uniform wind with the coupling flag unset, "projection", and "wrf", and checks: the default reproduces "projection" bit for bit, "wrf" tracks Rothermel's head rate to within 3%, and the default/"projection" head rate stays degraded toward the Wulff-shape tip -- reproducing the exact numbers validated for this PR (head_x(2100s) 3217.5 m default/projection vs 4091.6 m wrf). Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com>
…parse time, add FireDirectionalShape wrf deck, docs - Revert an accidental, undocumented change to the front-normal stencil in fill_directional_ros (one-sided/WENO derivatives instead of central differences) that slipped in during the earlier rebase and broke the "bit-for-bit identical to prior code" claim for the default coupling mode. Source/Fire/ERF_DirectionalRos.H now matches origin/ERF-Hazard exactly except for the intended wind_coupling additions. - Abort at parse time if directional_wind_coupling = "wrf" is combined with directional_shape = "ellipse" or a non-Rothermel ros_model (including hybrid), rather than silently no-oping. - Add a wrf deck and check to Exec/RegTests/FireDirectionalShape (point ignition), alongside the existing FireWrfWindCoupling (finite line). - Re-verified against the rebuilt binary: FireDirectionalShape 34/34 checks pass across 6 decks, with isotropic/projection/ellipse/ellipse_anderson reproducing their previously-recorded reference rates exactly. FireWrfWindCoupling 8/8 checks pass; with the corrected stencil the projection formula's head settles at 71.1% of Rf (not the 48.6% measured before the revert), wrf still reaches 100.0%, and default reproduces projection bit for bit (max|delta phi| = 0 over 71 shared plotfiles). - Document the four new keys in Docs/sphinx_doc/Inputs.rst, theory/fire_propagation.rst and inputs_fire_master_reference; add rothermel_ros_cell's missing @PARAM entries; position directional_wind_coupling = "wrf" against directional_shape = "ellipse" in the file comments; fix "circle" -> "stadium" support-function wording; note prescribed_wind is exact only with use_waf=false and use_terrain_wind=false. Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com>
Renamed so the option's name describes what it does rather than naming it after WRF-Fire specifically: since n |grad(phi)| = grad(phi), this coupling's R(n) turns the level-set equation into pure advection by a wind-aligned velocity plus isotropic growth -- the same reduction used in the original level-set fire-spread formulation (Mandel, Beezley & Kochanski 2011), which is where "advective" comes from. WRF-Fire's fire_ros happens to use the same functional form, and that cross-reference is kept in the comments/docs. - DIRECTIONAL_WIND_COUPLING_WRF -> _ADVECTIVE; erf.fire.directional_wind_coupling = "wrf" -> "advective" (ERF_DirectionalRos.H, ERF_FireParams.H, ERF_Rothermel.H) - Exec/RegTests/FireWrfWindCoupling -> FireAdvectiveWindCoupling, including its inputs_wrf -> inputs_advective deck, check/run scripts, and README - Exec/RegTests/FireDirectionalShape's wrf deck -> advective (inputs_wrf -> inputs_advective, check_firedirectionalshape.py's "wrf" kind -> "adv") - Docs/sphinx_doc/Inputs.rst, theory/fire_propagation.rst and inputs_fire_master_reference updated to match Re-verified against a clean rebuild: FireDirectionalShape 34/34 checks and FireAdvectiveWindCoupling 8/8 checks pass, with identical numbers to before the rename (e.g. advective head 0.24849 m/s / -0.38%, same as the prior wrf deck) -- confirms this is a pure rename with no behavior change. Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com>
|
Superseded by #421 — the branch this was on (ERF-Hazard) was renamed/moved out from under it, so rather than rebasing 20+ commits of drift I redid the same fix (plus everything from the review above) cleanly on top of the current ERF-Fire tip. Closing this one. |
2 participants
Summary
The current directional projection formula evaluates Rothermel's wind
factor with the wind projected onto the front normal before raising it
to the B power, R(theta) ~ (U cos theta)^B. This is non-convex once
phi_w(B-1) > 1, so a finite ignition line's endpoints grow a Wulff wedge
whose tip falls well below the head rate over time - verified against a
Munoz-Esparza-style finite-line test case, where the head rate degraded
from ~100% to ~40% of the theoretical Rf over 2100s.
erf.fire.directional_wind_coupling = "wrf" (default: "projection",
unchanged) instead exponentiates Rothermel's wind/slope factor from the
raw, unprojected wind speed and scales the result by the cosine to the
spread direction afterward, matching WRF-Fire's fire_ros
(module_fr_fire_phys.F). This keeps R(theta) linear in cos(theta) --
the convex support function of a circle -- so the wedge does not form;
the same test case now tracks the theoretical head rate to within 0.5%
for the full run.
Default behavior is verified bit-for-bit identical to the prior code
(max|delta phi| = 0 over the whole domain on the same test case).
Also includes two small pieces of test infrastructure used to validate
this:
interpolation with an exact constant wind, for controlled
idealized-wind test cases.
(default on, unchanged) -- needed when a test's wind speed exceeds the
cap.
Test plan
"wrf" tracks theoretical Rf to within 0.5% for the full 2100s run
(vs. degrading to ~40% of Rf with the default "projection" mode)
pre-change code